Dean Adams, Iowa State University
26 May, 2020
Biologists observe wondrous diversity in taxa and traits
Morphometricians wish to explain the evolution of phenotypic diversity
Our concern: how do we characterize patterns, and hypothesize processes?
Trait correlations often used to study coevolution and adaptation
Species values commonly utilized
The problem?
Species are not independent of one another
Phylogenetic comparative methods condition the data on the phylogeny during the analysis
Allows one to assess trait covariation while accounting for the non-independence due to shared evolutionary history
The conceptual development of PCMs
Phylogenetic regression/anova/association models
Phylogenetic Generalized Least Squares methods
Phylogenetic PLS
Phylogenetic ordination
Exploring and modelling evolutionary processes
Phylogenetic signal
Evolutionary rates
Evolutionary models
PCMs condition the data on the phylogeny under an evolutionary model
Data conditioning: Account for phylogenetic non-independence
Evolutionary model: How trait variation is expected to accumulate
Most PCMs use GLS (generalized least squares) as a model:
\[\small\mathbf{Y}=\mathbf{X{\hat{\beta}}+\epsilon}\] Here, \(\small\epsilon\) is not iid but is \(\small\sim\mathcal{N}(0,\textbf{V})\): containing expected covariation between taxa as described by the phylogeny (in matrix \(\small\mathbf{V}\))
\(\small\mathbf{V}\) is the phylogenetic covariance matrix
Describes the amount of evolutionary time species share via common ancestry (and thus how similar they are expected to be)
\[\small\mathbf{Y}=\mathbf{X{\hat{\beta}}+\epsilon}\]
Model design (\(\small\mathbf{X}\)) describes the type of analysis
Parameters of model (and model significance) obtained in various ways
Evaluate \(\small\mathbf{Y}=\mathbf{X}\mathbf{\beta } +\mathbf{E}\) in a phylogenetic context
The workhorse of PCMs
ANOVA and regression models that account for phylogeny
Requires a model of evolutionary change: typically Brownian motion (BM)
ANOVA and regression models that account for phylogeny
Requires a model of evolutionary change: typically Brownian motion (BM)
\(\Delta\mu = 0\)
\(\sigma^2\) (variance among taxa) \(\uparrow\) \(\propto\) time
ANOVA and regression models that account for phylogeny
Requires a model of evolutionary change: typically Brownian motion (BM)
\(\Delta\mu = 0\)
\(\sigma^2\) (variance among taxa) \(\uparrow\) \(\propto\) time
Statistical model: \(\small\mathbf{Y}=\mathbf{X}\mathbf{\beta } +\mathbf{E}\) where \(\small\mathbf{E} \sim \mathcal{N}(0,\textbf{V})\)
Several implementations possible (yield identical results if implemented properly)
Estimate contrast scores between pairs of taxa (tips or nodes)
Use contrasts for analyses (OLS solution)
Coefficients found as: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}_{pic} \mathbf{X}_{pic}\right )^{-1}\left ( \mathbf{X}^{T}_{pic} \mathbf{Y}_{pic}\right )\)
What is the association between Y and X?
## Registered S3 method overwritten by 'geiger':
## method from
## unique.multiPhylo ape
What is the association between Y and X?
## Analysis of Variance Table
##
## Response: pic.y
## Df Sum Sq Mean Sq F value Pr(>F)
## pic.x 1 14.3519 14.3519 19.285 0.00461 **
## Residuals 6 4.4651 0.7442
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pic.x
## 0.8846971
Use GLS model with non-independent error structure
Statistical model: \(\small\mathbf{Y}=\mathbf{X}\mathbf{\beta } +\mathbf{E}\) where \(\small\mathbf{E} \sim \mathcal{N}(0,\textbf{V})\)
Use GLS model with non-independent error structure
Statistical model: \(\small\mathbf{Y}=\mathbf{X}\mathbf{\beta } +\mathbf{E}\) where \(\small\mathbf{E} \sim \mathcal{N}(0,\textbf{V})\)
\(\textbf{V}\) is the phylogenetic covariance matrix
Describes expected covariance among taxa due to shared evolutionary history (typically under BM)
Coefficients found as: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}^{-1}\mathbf{Y}\right )\)
## Denom. DF: 6
## numDF F-value p-value
## (Intercept) 1 12.87792 0.0115
## X 1 19.28544 0.0046
## (Intercept) X
## -2.3296557 0.8846971
Identical results to PICs!
To understand what PGLS is doing, consider the following methods for obtaining model parameters
OLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}\mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{Y}\right )\)
To understand what PGLS is doing, consider the following methods for obtaining model parameters
OLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}\mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{Y}\right )\)
Unweighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1}\mathbf{Y}\right )\) where
\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right) \]
To understand what PGLS is doing, consider the following methods for obtaining model parameters
OLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}\mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{Y}\right )\)
Unweighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1}\mathbf{Y}\right )\) where
\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right) \]
Weighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{phy}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{phy}^{-1}\mathbf{Y}\right )\) where
\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} (v_1+v_2) & v_{12} & 0 \\ v_{12} & (v_1+v_2) & 0 \\ 0 & 0 & 1 \end{array} \right) \]
In PGLS, the weights are the phylogenetic distances, which describe the phylogenetic non-independence
To understand what PGLS is doing, consider the following methods for obtaining model parameters
OLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T}\mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{Y}\right )\)
Unweighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{id}^{-1}\mathbf{Y}\right )\) where
\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right) \]
Weighted GLS model: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}_{phy}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}_{phy}^{-1}\mathbf{Y}\right )\) where
\[\tiny\mathbf{V}_{id}= \left( \begin{array}{ccc} (v_1+v_2) & v_{12} & 0 \\ v_{12} & (v_1+v_2) & 0 \\ 0 & 0 & 1 \end{array} \right) \]
In PGLS, the weights are the phylogenetic distances, which describe the phylogenetic non-independence
Utilize OLS transformation of GLS models
Phylogenetic transformation matrix: \(\small\mathbf{T}=\left ( \mathbf{U}\mathbf{W}^{-1/2} \mathbf{U}^{T}\right )^{-1}\)
U and W are eigenvectors and eigenvalues of V
Transformed data: \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\)
Transformed design matrix: \(\small\tilde{\mathbf{X}}=\mathbf{TX}\)
Coefficients solved as: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{\tilde{X}}^{T} \mathbf{\tilde{X}}\right )^{-1}\left ( \mathbf{\tilde{X}}^{T} \mathbf{\tilde{Y}}\right )\)
## Df SS MS Rsq F Z Pr(>F)
## X 1 14.3519 14.3519 0.76271 19.285 1.8228 0.002 **
## Residuals 6 4.4651 0.7442 0.23729
## Total 7 18.8170
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Intercept) X
## -2.3296557 0.8846971
Identical results to PICs & PGLS!
PIC, PGLS, and Phylogenetic transform yield identical parameter estimates
\[\tiny \hat{\mathbf{\beta }}=\left ( \mathbf{X}^{T} \mathbf{V}^{-1} \mathbf{X}\right )^{-1}\left ( \mathbf{X}^{T} \mathbf{V}^{-1}\mathbf{Y}\right )\]
\[\tiny =\left ( \mathbf{X}^{T}_{pic} \mathbf{X}_{pic}\right )^{-1}\left ( \mathbf{X}^{T}_{pic} \mathbf{Y}_{pic}\right )\]
\[\tiny =\left ( \mathbf{\tilde{X}}^{T} \mathbf{\tilde{X}}\right )^{-1}\left ( \mathbf{\tilde{X}}^{T} \mathbf{\tilde{Y}}\right )\]
Model significance typically accomplished using parametric procedures
F-ratios compared to parametric F-distribution
Optimize \(\log{L}\) for model
\[\tiny \log{L}=\log{ \left(\frac{exp{\left({-\frac{1}{2}{\left({\mathbf{Y}-E(\mathbf{Y})}\right)^{T}}} \mathbf{V}^{-1}\left({\mathbf{Y}-E(\mathbf{Y})}\right)\right)}} {\sqrt{\left({2\pi}\right)^{Np}\times\begin{vmatrix} \mathbf{V} \end{vmatrix}}}\right)}\]
PROBLEM: Parametric significance testing suffers from Rao’s paradox
Power reduces as p-dimensions increase
Another solution required
Transform data: \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\), \(\small\tilde{\mathbf{X}}_{F}=\mathbf{TX}_{F}\), \(\small\tilde{\mathbf{X}}_{R}=\mathbf{TX}_{R}\)
Obtain Parameter estimates: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{X_{F}}}\right )^{-1}\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{Y}}\right )\), and permute residuals (as described previously)
Transform data: \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\), \(\small\tilde{\mathbf{X}}_{F}=\mathbf{TX}_{F}\), \(\small\tilde{\mathbf{X}}_{R}=\mathbf{TX}_{R}\)
Obtain Parameter estimates: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{X_{F}}}\right )^{-1}\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{Y}}\right )\), and permute residuals (as described previously)
Appropriate for ANOVA, regression and more complex PGLS models with high-dimensional data (Adams & Collyer. Evol. 2018)
RULE: Phylo-transform first, permute second!
Transform data: \(\small\tilde{\mathbf{Y}}=\mathbf{TY}\), \(\small\tilde{\mathbf{X}}_{F}=\mathbf{TX}_{F}\), \(\small\tilde{\mathbf{X}}_{R}=\mathbf{TX}_{R}\)
Obtain Parameter estimates: \(\small\hat{\mathbf{\beta }}=\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{X_{F}}}\right )^{-1}\left ( \mathbf{\tilde{X_{F}}}^{T} \mathbf{\tilde{Y}}\right )\), and permute residuals (as described previously)
Appropriate for ANOVA, regression and more complex PGLS models with high-dimensional data (Adams & Collyer. Evol. 2018)
RULE: Phylo-transform first, permute second!
Does head shape covary with body size among Plethodon salamander species?
Are body dimensions associated with overall size across Plethodon salamander species?
## Df SS MS Rsq F Z Pr(>Cohen's f-squared)
## svl 1 0.0006586 0.00065861 0.07039 3.0289 2.0186 0.012
## Residuals 40 0.0086977 0.00021744 0.92961
## Total 41 0.0093563
##
## svl *
## Residuals
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
YES! There is a significant association between head shape and body size in salamanders
Does head shape differ between high-elevation and low-elevation salamander species?
Does head shape differ between high-elevation and low-elevation salamander species?
## Df SS MS Rsq F Z Pr(>Cohen's f-squared)
## elev 1 0.0006142 0.00061420 0.06565 2.8103 1.8621 0.033
## Residuals 40 0.0087421 0.00021855 0.93435
## Total 41 0.0093563
##
## elev *
## Residuals
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
YES! Head shape differs in high- and low-elevation species
One can visualize group dispersion in morphospace…
#{r echo = FALSE, eval=TRUE,out.width="80%"} #plot.res <- gm.prcomp(shape,phy=plethtree, data=gdf) #ord.plot <- plot(plot.res,phylo = FALSE, pch=21, bg=gdf$elev, cex=1.5) #shapeHulls(ord.plot, groups = gdf$elev, # group.cols = c("red", "black"), # group.lwd = rep(1, 2), group.lty = c(2, 1)) #legend("topright", levels(gdf$elev), # col = c("black", "red"), # lwd = rep(1,2), lty = c(2, 1)) #
… and via TPS deformation grids
How groups are distributed on the phylogeny can make a difference: Adams and Collyer. Evol. (2018)
Too few group ‘state’ changes across phylogeny results in lower power